
// MPI + PETSc includes
#include "libmesh_project_soln.h"


// Global collections of MeshFunctions are necessary to work with
// global functions fptr and gptr.
// TODO: libMesh needs functor-based alternatives to these types of
// function arguments

/*
// Return the function value on the old mesh and solution
Number LibMeshProjection::fptr(const Point& p,
            const Parameters&,
            const std::string& libmesh_dbg_var(sys_name),
            const std::string& unknown_name)
{
	libmesh_assert(sys_name == current_sys_name);
	libmesh_assert(mesh_functions.count(unknown_name));
	libmesh_assert(mesh_functions[unknown_name]);

	this->projection_log->push("compute_value()");
	MeshFunction &meshfunc = *mesh_functions[unknown_name];
	Number value = meshfunc(p) ;
	this->projection_log->pop("compute_value()");

	return value;
}

// Return the function gradient on the old mesh and solution
Gradient LibMeshProjection::gptr(const Point& p,
              const Parameters&,
              const std::string& libmesh_dbg_var(sys_name),
              const std::string& unknown_name)
{
	libmesh_assert(sys_name == current_sys_name);
	libmesh_assert(mesh_functions.count(unknown_name));
	libmesh_assert(mesh_functions[unknown_name]);

	this->projection_log->push("compute_derivative()");
	MeshFunction &meshfunc = *mesh_functions[unknown_name];
	Gradient grad = meshfunc.gradient(p) ;
	this->projection_log->push("compute_derivative()");

	return grad;
}
*/

int LibMeshProjection::apply(std::string output_file)
{
	PetscPrintf(PETSC_COMM_WORLD, "\n[Starting the libMesh projection via MeshFunction]\n" );

	unsigned int n_systems = source->n_systems();
	libmesh_assert(target->n_systems() == n_systems);

	// For each system, serialize the solution so we can project it onto
	// a potentially-very-different partitioning
	for (unsigned int i = 0; i < 1/*n_systems*/; ++i)
    {
		System &old_sys = source->get_system(i);
		std::string current_sys_name = old_sys.name();

		std::cout << "Projecting system " << current_sys_name << std::endl;

		libmesh_assert (target->has_system(current_sys_name));

		System &new_sys = target->get_system(current_sys_name);
		unsigned int n_vars = old_sys.n_vars();
		libmesh_assert(new_sys.n_vars() == n_vars);

		this->projection_log->push("vectors_create()");

		AutoPtr<NumericVector<Number> > comparison_soln =
		NumericVector<Number>::build();
		std::vector<Number> global_soln;
		old_sys.update_global_solution(global_soln);
		comparison_soln->init(old_sys.solution->size(), true, SERIAL);
		(*comparison_soln) = global_soln;

		this->projection_log->pop("vectors_create()");

		this->projection_log->push("meshFunction_create()");

		// For each variable, construct a MeshFunction returning that
		// variable's value
		for (unsigned int j = 0; j != n_vars; ++j)
		{
			std::cout << " with variable " << old_sys.variable_name(j) << std::endl;

			MeshFunction *mesh_func =
			new MeshFunction(*source, *comparison_soln,
							 old_sys.get_dof_map(), j);
			mesh_func->init();
			mesh_functions[old_sys.variable_name(j)] = mesh_func;
		}

		this->projection_log->pop("meshFunction_create()");

		this->projection_log->push("project_solution()");

		// Project all variables to the new system
		new_sys.project_solution(new LibMeshProjection::FunctionValue(*this, current_sys_name));

		this->projection_log->pop("project_solution()");

		this->projection_log->push("meshFunction_destroy()");

		// Clean up the MeshFunctions here so we don't bloat memory
		for (unsigned int j = 0; j != n_vars; ++j)
		delete mesh_functions[old_sys.variable_name(j)];

		this->projection_log->pop("meshFunction_destroy()");

	}

  return 0 ;

}

